******************************************************************************************
* ONLINE APPENDIX FIGURE A6
*
* 1) ESTIMATE THE EFFECT OF THE VENEZUELAN IMMIGRATION ON WAGES, EMPLOYMENT, 
* UNEMPLOYMENT, AND LABOR FORCE PARTIPATION ON NATIVES BY SKILL LEVEL
*  
* LAST MODIFIED: JULY, 2025
******************************************************************************************

clear all 

* INSTALL SCHEMES FOR GRAPHS
*ssc install blindschemes, replace all 
*ssc install schemepack, replace /*based on https://github.com/asjadnaqvi/Stata-schemes*/

**# DIRECTORIES

* PAPER DIRECTORY DROPBOX
global fiscal "/Users/andres/Library/CloudStorage/Dropbox/2 Papers/2023/Paper_Fiscal_Costs"

* SUBDIRECTORIES
global data "$fiscal/1 Data/12 Final Datasets"
global programs "$fiscal/2 Programs/2 Figures and Tables paper"
global results "$fiscal/3 Results"

**# LOAD DATASET AND CREATE NEW VARIABLES

* LOAD DATASET WITH VARIABLES NEEDED
use "$data/LaborEffects.dta", clear

**#SAMPLE SELECTION 

* ONLY METROPOLITAN AREAS 
drop if area >= 81 /* not observations for these areas*/
drop if area == 25 /* not observations */
keep if area !=. 

* DUMMIES FOR YEAR AND MSA
tab year, gen(year)
tab area, gen(area)

* KEEP ONLY WORKING-AGE POPULATION
drop if (age<15|age>64) 

* KEEP IF NOT ENROLLED IN SCHOOL
*keep if att_school == 2 

* KEEP THOSE THAT WORK LESS THAN 100 HOURS A WEEK  
*drop if weekly_hours>100

* DELIMIT THE ANALYSIS FOR THE YEARS 2013-2018. 
drop if year > 2018

* DISPLAY FORMAT FOR SAMPLING WEIGHTS 
format %12.0f fex12

* GENERATE ROUNDED WEIGHTS FOR SPECIFIC COMMANDS THAT REQUIRE THEM 
gen fex12_round = fex12
replace fex12_round = round(fex12_round) 

**# GLOBAL MACROS FOR VARIABLES
global controls year2-year6 area2-area23 i.sex1 age c.age#c.age educ2-educ5
global x        mig_share2013
global z1       IV2005
global z2       IVdist
global w        fex12

**# ESTIMATION RESULTS

**## 1) EFFECT OF THE VENEZUELAN IMMIGRATION ON NATIVES' UNEMPLOYMENT

levelsof skill_group, local(skill)
foreach j of local skill {
	
display "This is skill `j'"

matrix A`j' = J(3,3,.) /* create matrix to save coefficients and confidence intervals*/
matrix colnames A`j' = b ll ul 
matrix rownames A`j' = OLS IV1 IV2 

**### OLS
ivreg2 unemp $x $controls if labforce == 1 & area != . & immig == 0 & skill_group == `j' [pw = $w], cluster(area)  partial(area2-area23) small 
matrix A`j'[1,1] = _b[$x]
matrix A`j'[1,2] = _b[$x] - invttail(e(df_r),0.025)*_se[$x]
matrix A`j'[1,3] = _b[$x] + invttail(e(df_r),0.025)*_se[$x]

**### IV1 (ENCLAVE INSTRUMENT)
ivreg2 unemp ($x = $z1) $controls if labforce == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix A`j'[2,1] = _b[$x]
matrix A`j'[2,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix A`j'[2,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]
**#### ANDERSON-RUBIN CI

**### IV2 (DISTANCE INSTRUMENT)
ivreg2 unemp ($x = $z2) $controls if labforce == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix A`j'[3,1] = _b[$x]
matrix A`j'[3,2] = _b[$x] - invttail(e(df_r),0.025)*_se[$x]
matrix A`j'[3,3] = _b[$x] + invttail(e(df_r),0.025)*_se[$x]

local ++j
}

matrix list A1
matrix list A2
matrix list A3


**## 2) EFFECT OF THE VENEZUELAN IMMIGRATION ON NATIVES' EMPLOYMENT

levelsof skill_group, local(skill)
foreach j of local skill {
	
display "This is skill `j'"

matrix B`j' = J(3,3,.)
matrix colnames B`j' = b ll ul 
matrix rownames B`j' = OLS IV1 IV2 

**### OLS
ivreg2 emp $x $controls if working_age == 1 & area != . & immig == 0 & skill_group == `j' [pw = $w], cluster(area)  partial(area2-area23) small 
matrix B`j'[1,1] = _b[$x]
matrix B`j'[1,2] = _b[$x] - invttail(e(df_r),0.025)*_se[$x]
matrix B`j'[1,3] = _b[$x] + invttail(e(df_r),0.025)*_se[$x]

**### IV1 (ENCLAVE INSTRUMENT)
ivreg2 emp ($x = $z1) $controls if working_age == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix B`j'[2,1] = _b[$x]
matrix B`j'[2,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix B`j'[2,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

**### IV2 (DISTANCE INSTRUMENT)
ivreg2 emp ($x = $z2) $controls if working_age == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix B`j'[3,1] = _b[$x]
matrix B`j'[3,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix B`j'[3,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

local ++j
}

matrix list B1
matrix list B2
matrix list B3


**##  3) EFFECT OF THE VENEZUELAN IMMIGRATION ON NATIVES' WAGES 
* NOTE: USUALLY THE SAMPLE IS RESTRICTED TO NATIVES IN MSAs WHO ARE 
* NOT ATTENDING SCHOOL AND NOT SELF-EMPLOYED. HOWEVER, MOST THE 
* PAPERS EXPLORING SIMILAR TOPICS FOR COLOMBIA DON'T DO SO.  

preserve
**### SAMPLE CONSIDERATIONS WAGE ESTIMATIONS
* KEEP IF NOT SELF-EMPLOYED 
*keep if p6430 != 4

* KEEP ONLY SALARIED WORKERS
*keep if p6430 == 1

* KEEP THOSE THAT WORK LESS THAN 100 HOURS A WEEK  
*drop if p6800>100

* HANDLING OUTLIERS
levelsof year, local(year)
foreach j of local year {
display "this is year `j'"
centile hourly_real_wages if year == `j', cent(0.5 99.5)
return list 
local cent1 = r(c_1)
display `cent1'
local cent99 = r(c_2)
display `cent99'
* keep only de wage distribution between 1 and 99
drop if hourly_real_wages > `cent99'  
drop if hourly_real_wages < `cent1' 
}

levelsof skill_group, local(skill)
foreach j of local skill {
	
display "This is skill `j'"

matrix C`j' = J(3,3,.)
matrix colnames C`j' = b ll ul 
matrix rownames C`j' = OLS IV1 IV2 

**### OLS
ivreg2 log_hourly_real_wages $x $controls if emp == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], cluster(area)  partial(area2-area23) small 
matrix C`j'[1,1] = _b[$x]
matrix C`j'[1,2] = _b[$x] - invttail(e(df_r),0.025)*_se[$x]
matrix C`j'[1,3] = _b[$x] + invttail(e(df_r),0.025)*_se[$x]

**### IV1 (ENCLAVE INSTRUMENT)
ivreg2 log_hourly_real_wages ($x = $z1) $controls if emp == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area)  partial(area2-area23) small 
matrix C`j'[2,1] = _b[$x]
matrix C`j'[2,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix C`j'[2,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

**### IV2 (DISTANCE INSTRUMENT)
ivreg2 log_hourly_real_wages ($x = $z2) $controls if emp == 1 & immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area)  partial(area2-area23) small 
matrix C`j'[3,1] = _b[$x]
matrix C`j'[3,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix C`j'[3,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

local ++j
}
restore /*restore sample before considering wage outliers*/

matrix list C1
matrix list C2
matrix list C3


**## 4: EFFECT OF THE VENEZUELAN IMMIGRATION ON NATIVES' LABOR FORCE PARTICIPATION

levelsof skill_group, local(skill)
foreach j of local skill {
	
display "This is skill `j'"

matrix D`j' = J(3,3,.)
matrix colnames D`j' = b ll ul 
matrix rownames D`j' = OLS IV1 IV2 

**### OLS
ivreg2 labforce $x $controls if immig == 0 & area != . & skill_group == `j' [pw = $w], cluster(area)  partial(area2-area23) small 
matrix D`j'[1,1] = _b[$x]
matrix D`j'[1,2] = _b[$x] - invttail(e(df_r),0.025)*_se[$x]
matrix D`j'[1,3] = _b[$x] + invttail(e(df_r),0.025)*_se[$x]

**### IV1 (ENCLAVE INSTRUMENT)
ivreg2 labforce ($x = $z1) $controls if immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix D`j'[2,1] = _b[$x]
matrix D`j'[2,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix D`j'[2,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

**### IV2 (DISTANCE INSTRUMENT)
ivreg2  labforce ($x = $z2) $controls if immig == 0 & area != . & skill_group == `j' [pw = $w], first cluster(area) partial(area2-area23) small 
matrix D`j'[3,1] = _b[$x]
matrix D`j'[3,2] = _b[$x ] - invttail(e(df_r),0.025)*_se[$x]
matrix D`j'[3,3] = _b[$x ] + invttail(e(df_r),0.025)*_se[$x]

local ++j
}

matrix list D1
matrix list D2
matrix list D3


**# PLOTTING EFFECTS

**## WAGES
*  T-BASED CIs
levelsof skill_group, local(skill)
foreach j of local skill {
coefplot (matrix(C`j'[,1]), ci((C`j'[,2] C`j'[,3])) ///
ciopts(recast(rcap))), ///
yline(0) vertical ///
coeflabel(1 r1 = "OLS" ///
		  2 r2 = "IV: enclave" ///
		  3 r3 = "IV: distance" , labsize(small)) ///
ytitle("Effect on log hourly wages skill `j'", size(medsmall))  ///
graphregion(color(white)) scheme(plottig) ///
		title("A. Wages - Skill `j'", size(medium))
		 tempfile resultsC`j'
		 graph save "`resultsC`j''"		
}	

		 		 
**## UNEMPLOYMENT 		 
* T-BASED CIs
levelsof skill_group, local(skill)
foreach j of local skill {
coefplot (matrix(A`j'[,1]), ci((A`j'[,2] A`j'[,3])) ///
ciopts(recast(rcap))), ///
yline(0) vertical ///
coeflabel(1 r1 = "OLS" ///
		  2 r2 = "IV: enclave" ///
		  3 r3 = "IV: distance" , labsize(small)) ///
ytitle("Effect on P(unemployment = 1) skill `j'", size(medsmall))  ///
graphregion(color(white)) scheme(plottig) ///
title("B. Unemployment skill `j'", size(medium))
tempfile resultsA`j'
graph save "`resultsA`j''"
}


**## EMPLOYMENT
* T-BASED CIs
levelsof skill_group, local(skill)
foreach j of local skill {
coefplot (matrix(B`j'[,1]), ci((B`j'[,2] B`j'[,3])) ///
ciopts(recast(rcap))), ///
yline(0) vertical ///
coeflabel(1 r1 = "OLS" ///
		  2 r2 = "IV: enclave" ///
		  3 r3 = "IV: distance" , labsize(small)) ///
ytitle("Effect on P(employment = 1) skill `j'", size(medsmall))  ///
graphregion(color(white)) scheme(plottig) ///
		title("C. Employment skill `j'", size(medium))
		 tempfile resultsB`j'
		 graph save "`resultsB`j''"		 
}

		 
**## LABOR FORCE PARTICIPATION
* T-BASED CIs
levelsof skill_group, local(skill)
foreach j of local skill {
coefplot (matrix(D`j'[,1]), ci((D`j'[,2] D`j'[,3])) ///
ciopts(recast(rcap))), ///
yline(0) vertical ///
coeflabel(1 r1 = "OLS" ///
		  2 r2 = "IV: enclave" ///
		  3 r3 = "IV: distance" , labsize(small)) ///
ytitle("Effect on P(labor force = 1) skill `j'", size(medsmall))  ///
graphregion(color(white)) scheme(plottig) ///
		title("D. Labor force skill `j'", size(medium))
		 tempfile resultsD`j'
		 graph save "`resultsD`j''"
}


* (T-BASED CIs) -> COMBINE GRAPS PER OUTCOME TO FACILITATE Y AXIS SCALE INTERPRETATION 
* WAGES
graph combine "`resultsC1'" "`resultsC2'" "`resultsC3'" ///
, graphregion(color(white)) cols(3) rows(1) xcommon ycommon imargin(0 0 0 0)  ///
 iscale(0.72) l2title("Effect on log hourly wages") scheme(plottig) 
 tempfile wages
 graph save "`wages'"

* UNEMPLOYMENT  
graph combine  "`resultsA1'" "`resultsA2'" "`resultsA3'" ///
, graphregion(color(white)) cols(3) rows(1) xcommon ycommon imargin(0 0 0 0)  ///
 iscale(0.72) l2title("Effect on P(unemployment = 1)") scheme(plottig) 
 tempfile unemployment
 graph save "`unemployment'" 

* EMPLOYMENT
graph combine "`resultsB1'" "`resultsB2'" "`resultsB3'" ///
, graphregion(color(white)) cols(3) rows(1) xcommon ycommon imargin(0 0 0 0)  ///
 iscale(0.72) l2title("Effect on P(employment = 1)") scheme(plottig) 
 tempfile employment
 graph save "`employment'" 
 
* LABOR FORCE PARTICIPATION 
graph combine "`resultsD1'" "`resultsD2'" "`resultsD3'" ///
, graphregion(color(white)) cols(3) rows(1) xcommon ycommon imargin(0 0 0 0)  ///
 iscale(0.72) l2title("Effect on P(labor force = 1)") scheme(plottig) 
  tempfile laborforce
 graph save "`laborforce'" 
  
**## FIGURE FOR THE ONLINE APPENDIX: LABOR MARKET EFFECTS BY SKILL (T-BASED CIs) 
graph combine  "`wages'" "`unemployment'" "`employment'" "`laborforce'", graphregion(color(white)) cols(1) rows(4) xcommon ycommon imargin(0 0 0 0) play(editOutcomesBySkill) ///
iscale(0.72) scheme(plottig)  
qui graph save "$results/1 Figures Paper/FigA6_1.gph", replace 
qui graph export "$results/1 Figures Paper/FigA6_1.eps", as(eps) replace		 
 

